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ABSTRACT 

We present the secular theory of coplanar -planet system, in the absence of mean motion 
resonances between the planets. This theory relies on the averaging of a perturbation to the 
two-body problem over the mean longitudes. We expand the perturbing Hamiltonian in Taylor 
series with respect to the ratios of semi-major axes which are considered as small parame- 
ters, without direct restrictions on the eccentricities. Next, we average out the resulting series 
term by term. This is possible thanks to a particular but in fact quite elementary choice of 
the integration variables. It makes it possible to avoid Fourier expansions of the perturbing 
Hamiltonian. We derive high order expansions of the averaged secular Hamiltonian (here, up 
to the order of 24) with respect to the semi-major axes ratio. The resulting secular theory is 
a generalization of the octupole theory. The analytical results are compared with the results 
of numerical (i.e., practically exact) averaging. We estimate the convergence radius of the 
derived expansions, and we propose a further improvement of the algorithm. As a particular 
application of the method, we consider the secular dynamics of three-planet coplanar system. 
We focus on stationary solutions in the HD 37124 planetary system. 

Key words: celestial mechanics - secular dynamics - analytical methods - stationary solu- 
tions - extrasolar planetary systems - stars: HD 37124 



1 INTRODUCTION 

The recent discoveries of extrasolar planetary systems bring new 
and interesting problems regarding their dynamical stability and 
long-term evolution. At present, at least 30 multi-planet systems 
are known and their number is still growing, thanks to refined tech- 
niques of observations. Surprisingly, the orbital parameters of these 
systems are very different from those typical in the Solar System 
architecture — large planetary masses and eccentricities are com- 
mon. Simultaneously, these systems usually are compact. Likely, 
this property is a consequence of the observational selection. The 
most effective detection techniques, like the radial velocity ob- 
servations, rely on indirect effects of mutual interactions between 
planets and their host star. Many of the known multi-planet sys- 
tems are supposed to be involved in short-term mean motion res- 
onances (MMRs). However, there are also configurations with rel- 
atively well separated orbits. In that case the secular interactions 
may lead to interesting dynamical phenomena. 

To study the long term dynamics of planetary systems, dif- 
ferent analytical and numerical techniques are used. The analytical 
approach is much more effective in the investigations of global, 
qualitative dynamics than widely applied numerical techniques (in- 
cluding fast indicators or direct numerical integration of the equa- 
tions of motion). The numerical experiments provide only limited 
(or local) information on the dynamical features of the studied con- 
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figurations. Usually, the interpretation of the results of massive cal- 
culations can be problematic without solid theoretical background. 
In contrast, analytical techniques offer much deeper insight into 
qualitative properties of motion. The analytical approach makes it 
possible to explore large volume of the phase space. This is crucial 
for the dynamical studies of extrasolar planetary systems detected 
during short time of observations. These observations have rela- 
tively large errors, they are typically irregularly sampled or degen- 
erated (in the sense that they can provide only limited information 
on the system state, like the radial velocity technique). This leads 
to poorly determined or unconstrained orbital and physical param- 
eters of the detected systems. In that case, the analytical theories 
help us to investigate and/or to detect global properties of the solu- 
tions to the equations of motion in observationally permitted ranges 
of the parameters. We can investigate in detail certain families of 
these solutions, their bifurcations and stability. Examples of such 
solutions are stationary solutions (equilibria) or periodic orbits that 
build a skeleton of the phase space. Investigating these families, we 
follow the classic methodology invented by Poincare. The analyti- 
cal theories are milestones for detailed numerical studies of partic- 
ular aspects of the dynamics. Hence, their constant development is 
always desirable. 

Moreover, due to extreme parameters of the studied con- 
figurations, the classic planetary theory developed so far is of- 
ten too week. For instanc e, the classic Lagrange-Laplace theory 
dMurrav & Dermottll2000h designed as a model of the secular dy- 
namics of planets in the Solar System fails in the case of large 
eccentricities and inclinations. Hence, new analytical and semi- 
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analytical theories are recently developed, breaking the limitations 
of the classic approach. One of the most effective techniques for 
studying the secu lar dynamics of extrasolar syste ms has been re- 
cently invented by Michtchenko & MalhotraU2004h and further de- 
veloped in ( Mic htchenko et alj|2006h . These papers are devoted to 
a study of two-planet configurations. In this work, we consider an 
analytical secular theory of a coplanar system of N-planets (point 
masses) under assumption that the orbital configurations are not 
involved in strong mean-motion resonances and that they are far 
from collisions zones of orbits. We calculated the averaged per- 
turbation in the form of power series with respect to the semi- 
major axes ratios up to very high order (equal to 24 in the present 
work). These expansions have no explicit limits on the eccentric- 
ities provided that the non-resonance condition is satisfied. Our 
development is elementary and is based on very basic properties 
of the Keplerian motion. Although it concerns the two-planet sys- 
tem, we show that it can be easily generalized for the case of N- 
planet configurations. Hence, the theory can be regarded not only 
as an attemp t to improve the secular theories for two-planet sys- 
tems in (e.gjRodnguez & Gallardol20 05 ; Henra rd~& Libertl 2005; 
Libert & Henrardl I2005L 120061 : |ji et all 120071 : IVeras & Arm itage 



20071) . relying on the cl assic expansion of the perturbing Hamil- 



tonian in eccentricities (Mur ray & Dermoid [2000; Ellis & Murrav 
2000) or the octupole theory that makes use on the averaging of 
the low-order ex pansion of the perturbing function in the semi- 
major axes ratio jFord et alJkood : iBlaes et alj|2002l ; iLee & Pealel 
120031) . We try to reduce the limitations of the classic theory of non- 
resonant systems. 

The plan of this paper is as follows. In Section 2, we describe 
a general model of a coplanar configuration of N planets. We in- 
troduce the expansion of perturbing Hamiltonian and a very sim- 
ple and basic algorithm of its averaging. We compare the results 
of the method with the outcome of the octupole theory, and we 
discuss some subtle differences between these theories. We also 
present the results of the tests of the expansion, taking as exam- 
ples a few known multi-planet configurations which apparently fit 
well in the framework of the non-resonant secular theory. The ex- 
act semi-numerical method is helpful to determine absolute bounds 
of the validity of the analytic approach. We also outline a further 
improvement of the averaging algorithm. In Section 3 we construct 
the secular model of the three-planet system and we perform a pre- 
liminary study of the secular dynamics of a few known three-planet 
configurations. In particular, we focus on the equilibria in the sec- 
ular problem, and we found i nteresting stationary solutions in the 
extras olar system HD 37124 dVogt et aT]|2005l ; [Gozdziews ki et al.l 
l2008h . 



2 THE SECULAR DYNAMICS OF A MULTI-PLANET 
SYSTEM 

The Hamiltonian of a multi-planet system with respect to 
canonical Poincare variables ( see, e.g.. lLaskar & Robutellll995l ; 
Michtchenko & Mal hotrj|200l) can be expressed by a sum of two 
terms, 



& — ^4epl + ^pcrti 



where 



£ V 2\h n 



(l) 



(2) 



is for the integrable part comprising of the direct sum of the rel- 
ative, Keplerian motions of N planets and the host star. Here, the 
dominant point mass of the star is mo, and m; <C mo, i = l,...,N 
are the point masses of the Af-planets. For each planet-star pair we 
define the mass parameter /j,- = k 2 (mo + mi) where k is the Gauss 
gravitational constant, and (3/ = (1/m; + 1 /mo) -1 are the so called 
reduced masses. Due to mutual interactions between the planets, 
the Keplerian part is perturbed by a function J{ pETt , 



N-l N 
i=l ;>! 



k 2 mimj p, Pj 
Ay mo 

direct part indirect part 



(3) 



where r,- are for the position vectors of the planets relative to the 
star, p,- are for their conjugate momenta relative to the barycenter of 
the whole (N+ l)-body system, Ay = ||r,- — rj\\ denote the relative 
distance between planets ( and j. It is well known that even in the 
simplest case of three point masses (the star and two mutually inter- 
acting planets), the problem is non-integrable and is not possible to 
obtain its exact analytical solutions. In practice, such solutions can 
be only derived in the form of approxi mations derived by means 
of different perturbation techniques ( e.g., lMurrav & Dermottl2000l ; 
lMorbidellill2002l; lFerraz-Melloll2007l) . 

To apply the canonical perturbation theory, we first transform 
9l to the following form: 



#(1,0) =^4epl (I) + ^per, (1,0), 



(4) 



where (1,0) stand for the action-angle variables, and #" pert (I,0) ~ 
£^4 cp i(I), where e <C 1 is a small parameter. In the absence of this 
perturbation, the system is trivially integrable. However, with the 
perturbation added, the dynamics of the full system become ex- 
tremely complex. In the realm of the Hamiltonian canonical the- 
ory, the approximate, analytical solutions to this problem may be 
derived by an expansion of the perturbation with respect to the 
small parameter and by subsequent simplification of the lowest or- 
der terms by means of appropriate canonical contact transforma- 
tions. This idea of Delaunay appears in many "incarnations". One 
of its first novel reali zations is known as the von Zeipel method 
(e.g.. lBrumberg| [l995). Much more improved version of this tech- 
nique that not require se ries inversion has been invented by iHoril 
( 1966) and next refined by |Depritl ( fl969l) . This theory is well known 
in the literature as the Lie-Hori-Deprit metho d. For an exce l lent re - 
view of these methods see a monograph by iFerraz-MeH ol d2007l) . 
Another method of seeking for approximate solutions to Eq.[4]re- 
lies d irectly on the averaging proposition (see, e.g., Arn old et al.l 
Il993h . Usually, the canonical angles can be divided onto two 
classes: fast and slow ones. By averaging the perturbing part with 
respect to the fast angles over their periods, we obtain the secular 
perturbation Hamiltonian which does not depend on these fast an- 
gles. Simultaneously, their conjugate momenta become integrals of 
the secular problem. In the planetary model with a dominant stellar 
mass, we have two natural time-scales of motion: the orbital mo- 
tion of the planets and a slow evolution of their orbits. Assuming 
that no strong mean motion resonances are present, and the system 
is far enough from collisions, the averaging makes it possible to re- 
duce the number of the degrees of freedom and to obtain qualitative 
information on the long-term changes of the slowly varying orbital 
elements (i.e., on the slow angles and their conjugate momenta). 

To apply each one of these methods, the Hamiltonian of 
the iV-planet system should be first transformed to the required 
form, Eq. [4] This can be accomplished by expressing it with re- 
spect to the canonical Poincare elements dMurrav & Dermot feOQOl : 
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Mich tchenko & Malhotr j| 2004b : 



// = A,;, 
gi = — OJi, 
= -Sli, 




Gi=mi-Jl-ej), 



(5) 



Hi=LiJl-ef (l-cos/j) 



where A,/ are the mean longitudes, a, stand for canonical semi-major 
axes, e,- are for the eccentricities, denote inclinations, 03/ are the 
longitudes of pericenter, and £2,- denote the longitudes of the as- 
cending node. We note that the transformation between the canon- 
ical orbital elements of Poincare, a,, e,-, /,-, B3;, £2; and associated 
Cartesian coordinates and momenta may be derived by the formal 
two-body transformation between classic (astro- centric) Keplerian 
elements and the Cartesian coordinates (e.g., iMorbidellj 120021 ; 
iFerraz-Mello et alii 20061) . Moreover, in the settings adopted here, 
the rectangular coordinates and momenta are understood through 
the Cartesian positions of planets relative to the star, and, according 
to the definition of the Poincare variables, the canonical momenta 
are taken relative to the barycenter of the system. 

The W-body Hamiltonian expressed in terms of the Poincare 
variables has the form of: 



•H 



'h 24 



■ ?(pea{Li,k,Gi,gi,Hi,hi) 



!=!,.. .JV 



In this Hamiltonian, // play the role of the fast angles. In the absence 
of strong MMRs, these angles can be eliminated by the following 
averaging formulae: 



1 



(2n) N 



i-2% r2% 
/ •■•/ ' 

Jo Jo 



(6) 



(As we see below, it can be also applied for some selected pairs 
of planets). Hence, the conjugate momenta L, become integrals of 
the secular model and the Keplerian part is a constant that does not 
contribute to the equations of motion. However, the calculation of 
the multiple integral is quite a difficult task which is a central part 
of the problem. We try to solve it with quite basic mathematical 
properties of the Keplerian osculating orbits. 

We should keep in mind that the averaging of the secular 
Hamiltonian in the problem over Keplerian motions implies trun- 
cation of the perturbation to first order in the masses (more gener- 
ally, to the first order in the perturbation parameter £). For large 
mutually interacting planets (or binary stars), the deviations of 
the true orbits from the Keplerian approximation during the or- 
bital period may become significant, and the secular theory may 
fail. Nevertheless, it is a common drawback of the idea behind 
Eq.|6] For the same reason, the classic perturbation techniques usu- 
ally fail in the case of close encounters between the planets. In 
such an instance, some other criteria helping to explore the sta- 
bility r egions may be applied, for example, the Hill stability crite- 
rion iMarchal & B ozis 1982; Gladman 1993; Ba rnes & Greenbergl 
l2006l : lMichtchenko et a l. 2008; Barnes & Greenberdl2008b~ 



2.1 The indirect part of the disturbing function 

We start with the averaging of the indirect part of the disturb- 
ing Hamiltonian. The resu lt of this averaging can be found in 
Brou wer" & Clemencd J 196 lh . nevertheless, to make this paper self- 
consistent, we present the calculations in detail. The indirect part 



is a scalar product of canonical momenta p,, which have the form 
of: 



, . , niimj . 



(7) 



where p" ; = [1 /m, + 1/(M — m,-)] ~ and M is the total mass of the 
system. The scalar product p,- ■ pj includes terms of the type of 
r, • tj. Moreover, each product p,- • p ; depends on all astro-centric 
velocities of the planets, r; (i = l,...,N). Apparently, to average 
out the indirect part of the disturbing function, we must compute 
multiple integral over all mean longitudes A,- (i = 1, . . . ,N). In fact, 
this integral can be reduced to a sum of double integrals computed 
for all pairs of planets, i.e., we average out expressions of the form 
of r,- • r ,-. The result is the following: 



1 



r2% i-2% 

I I r; • Ti dMid'M j = 5; ; 1 
JO Jo 



{2nf 



(8) 



where n,- denote mean motions of the planets and 8, j stands for the 
Kronecker delta. Note, that the averaging over the mean anoma- 
lies gives the same results as the averaging over the mean longi- 
tudes under the condition that the integration limits are set to and 
271, respectively. Clearly, the indirect part of the disturbing func- 
tion does not contribute to t he secular dynamics of the s ystem be- 



cause it depends on L,- only 1 Brouwer & Clemencd 196 lb . (see also 
iMichtchenko & Malhotr aj2004b . We note that this result is exact as 
far as the assumptions of the averaging principle are fulfilled (we 
are far enough from the MMRs and there are present two different 
time-scales in the problem). 



2.2 The direct part of the disturbing function 

Now we have a more difficult problem to resolve. For the iV-planet 
system, the averaged direct part of the disturbing function has the 
form of: 



N-\ N 
[=1 j>i 



(9) 



where the multiple integral Eq. l[6j over all mean anomalies is re- 
duced to a sum of secular Hamiltonians describing mutual interac- 
tions between all pairs of planets; i.e., for each pair where 
i < j and ci[ < aj, we have: 



(27t)- 



1 f 2% f 2K k 2 niimj 



lo Jo 



dMidMj. (10) 



Hence, the secular model of ,/V-planet system can be reduced to a 
simple sum of two-planet Hamiltonians. 

Now, we compute the double integral for a fixed pair of planets 
i and /. The distance between these planets is the following: 



A ''-> = y r f+ r j- inrjco&Vij, (ii) 

where yfu is the angle between vectors r; and rj: 
Ti-rj _ x/Xj+yiyj 



cos j — 

•i'l 

This formulae can be rewritten to: 
A; 



't'J 



. U = rjj 1-2- (x^+yA 



(12) 



(13) 



According to the Kepler problem theory, r; and rj may be expressed 
through the eccentric anomaly, E, or, equivalently, through the true 



© 2008 RAS, MNRAS OOO.QJfl]] 



4 C. Migaszewski and K. Gozdziewski 



anomaly, /. We write down appropriate expressions for planet (' and 
j, respectively: 



aj(l-e 2 j) 



+ ejcosfj 

Here, Ej is the eccentric anomaly of the inner planet in the selected 
pair of interacting bodies, and fj is the true anomaly of the outer 
planet. Moreover: 

— = cos (fj + AG5;j) , — sin (fj + ADJy) , 
where AGS,- = ((55; — GJ,-), and 



Xj = cij (cos Ej — e, ) , >',- = a, \J 1 — ef sin . 

The dependence of the two-body formulae on AG3,j may seem 
strange. However, when we investigate the co-planar system of two 
particular planets, we are free to choose the reference frame be- 
cause the mutual interaction between these planets depend only 
on their relative orbital phases. To be more specific, we must cal- 
culate the distance between planets i and j, or the scalar product 
r,- • r; = r,rj cos \|/,- t (Eq.ll It. That is obviously independent on the 
reference frame. In general, we could write: 



r,- • r = AjTj 



(14) 



where f/j are the orbital reference frames of the inner and outer 
planet, respectively, matrices A,-, Ay represent Eulerian rotations of 
flj to a common reference frame ( J ) for both orbits. Because this 
frame may be chosen freely, we fix the x-direction of the common 
frame along the apsidal line of the inner planet (still, for a particu- 
lar pair of planets). Then A,- = E and f j must be rotated by angle 
AG5; That can be repeated for each pair of planets in multi-planet 
system because the secular Hamiltonian is represented by a sum of 
formally independent two-planet terms. 

Finally, the inverted distance between planets i,j can be ex- 
pressed as follows: 

1 l+ejco&f] 
where a,j = at/ai < 1, and 



lAajj-lBOij + i 



-1/2 



(1 + ej cos//) 2 (l — e, cos£,' 
1 + ej cos fj 



1 



^(cos£,- — 6j) cos (fj + AB5; 



(15) 



(16) 



+ J 1 - ef sin E; sin (fj + AWij) 



Now we underline that the position of the inner planet is given 
through the eccentric anomaly while the position of the outer planet 
is given with respect to the true anomaly. The formulae under the 
square root are expressed by a polynomial of trigonometric func- 
tions and, as we can see below, that is critically important property 
making it possible to calculate the integral in Eq.llOl 

Now, we expand the inverse of the distance between planets 
i and j in Taylor series with respect to small parameter a, j. The 
series are evaluated around rx, j = as follows: 



1 



l + eycos/, 



(1- 



1 d'v 



II da 



where 



Aaf 



2Ba, 



hi ' 



a, :=o 



-1/2 



(17) 



(18) 



As the final result of this expansion, we obtain a polynomial with 
respect to trigonometric functions of the anomalies which has the 
general form of: 



1 



: £ [Cj,(cosEi)P< (sin£i) K (cos /,)« (sin fj) 



(19) 



Here, p = (pi,P2,P3,P4) £ N 4 is a vector of natural numbers 
and C p are coefficients depending on eccentricities and semi-major 
axes. One more step is still necessary. The integral in Eq.UOImust 
be computed with respect to the mean anomalies. Here, the clas- 
sic theories make use on sin/ and cos / expressed through Fourier 
series of the mean anomalies with coefficients dependent on the 
eccentricities. 

However, we found that it is possible to avoid these ex- 
pansions. Again, using the basic formulae of the Keplerian mo- 
tion, we perform a formal change of variables in Eq. [10] with: 

d Mi = IjdEj and d 9Aj = Jjdfj , where functions // = (Ei,ef) and 
3j = Jj(fj,ej) are defined with: 

3/2 



h(Ei,ei) = 1 - e,cos£ ; , fj(fj,ej) 



(20) 



(1 + e j cos fj) 

After this change of variables, the average in Eq. [TU]is equivalent 
to calculation of the following double integral: 
1 f2jl 



where A 



(2n) 



V 7 

■It) 2 Jo Jo 



271 2 1 

k minij — Ij (Ei ,et)j j (fj ,ej) dE t df j , 

(21) 



\j (ai,aj,ei,ej,Ei,fj). 



Fortunately, functions /; are again trigonometric polynomials 
and they do not change the general, polynomial form of Eq. [19] 
However, the second scaling function is not a polynomial with re- 
spect to cos/,- or sin/. Yet Eq.fTTlinvolves afactor (l +ej cos/y) . 
It cancels out one power of ( 1 + ej cos/;) appearing in the denom- 
inator of Eq. [20] In order to calculate the expansion in Eq. [T7] for 
(>0we differentiate D (a,- j) with respect to a,- j as the compos- 
ite function (Eq.ll8t. This operation emerges factors of the type of 
A r B s , where n = r + s ^ 1. Looking at the general form of A and 
B we see that the term (l +eycos/)) appears with natural powers 
larger than 1 and it cancels out remaining ( 1 + ej cos/y) in the de- 
nominator of Eq.[20] In this way, the general form of trigonometric 
polynomial in Eq.[l9]is preserved. Still, the free term in the Taylor 
expansion of © leads to an expression involving (l +ejCOsfj) '. 
Fortunately, we must integrate such term with limits from to 2% 
and this effectively can be reduced to averaging out rj over the or- 
bital period (or the whole range of the true anomaly). 

Finally, we can integrate Eq. [2T]term by term. Basically, the 
problem has been reduced to the calculation of definite integrals 
from products of trigonometric functions sin(x) and cos(x) in some 
natural powers. These integrals can be derived quite easily, at least 
in principle, nevertheless with increasing order of the Taylor ex- 
pansion, the calculations become extremely tedious. To accomplish 
them, we used MATHEMATIC A and fast AMD-Opteron computer. 

The final result of the averaging is an expansion of the secular, 
two-body Hamiltonian for a chosen pair of planets i and j: 



rtf (ij) = _ 



k 2 m,mj 



1 + a/1 



1=2 



l-ei 



^ iJ) (ei, 



,A05y) 



(22) 



Looking at the general form of this secular Hamiltonian, we learn 



© 2008 RAS, MNRAS 000.ITHT41 



Secular planetary problem 5 



that the role of a formal parameter in the power-series expansion of 
^/s K plays (apparently) the following expression: 



1-e 



ay(l 



•4) 



(23) 



Obviously, these series cannot converge if Xu ^ 1, hence we re- 
quire that Xjj < 1. Of course, this is only the necessary (and as we 
see below, very rough) condition for the convergence of these se- 
ries. However, as we explain in Sect. 2.5, attributing to X/j the role 
of a parameter deciding on the convergence of these series is in fact 
misleading because their divergence flows from quite a different 
source. 

A few first terms of the expansion of the secular Hamiltonian, 
Eq. [22] are listed below. The free term (for / = 0) is constant be- 
cause it depends on the mean semi-major axis aj only. The term 
with / = 1 vanishes identically. 

Terms of order 2 and 3 may be identified wi th the quadrupole 
and octupole secular Hamiltonian, respectively ( iFord et al.ll2000l ; 
iLee & Pealell2003h : 



3L-. 



(U) 



1 



3e- 



Higher order terms are the following: 



4 J) = 



1024 



^ u = _ 



[70 (ef +2) e?e]cos(2ACJy) 4 
ISef + 40ef + s) (3e 2 + 2)], 

[7 (3ef + 8) e)ef cos(3AC3, J ) 



105 
4096 I 



(24) 
(25) 

(26) 
(27) 



-2 (5 (ef + 4j ef + 8) lie] + 4j e,e ; cos(AG5,j 



65536 



[2079 (lef + 10) efejcos(4A0J, )7 -) + (28) 



+630 (l5ef + Wef + 48J [ej+2j e?^cos(2A05 i) y) 
+ 10 (35ef + 2\0ej + 168ef + 16) (l5e) + 40e 2 + 8 



We computed the expansion up to the order of 24. This expan- 
sion is available on the request in the form of raw MATHEMAT- 
ICA input file; also available in the form of on-line material after 
publishing this paper. 



2.3 A comparison with the octupole theory 

Here, we compare the results provided by the secular theory derived 
in the previous section with the results obtained wi th the help of the 
octupole theory of two planets dLee & Peale 2003) which has been 
obtained through averaging the perturbation Hamiltonian with the 
help of von Zeipel method up to the third order in 04 2 = K- It has 
been applied to study qualitative features of the secular dynamics in 
hierarchical planetary systems (i.e . with small a). A similar theory 
has been developed bv lFord et alj ( feOOOl) who investigated secular 
dynamics in hierarchical triple stellar systems with large separation 
of the third body. 

To make our d iscussion more tran sparent, we use, in this sec- 
tion, the notation of lLee & Pealel l l2003ri . Their equations (13), (14) 



and (15) have the following form: 



L 2 
Gi 



matrix 
mo +m\ 
(mo+mi)m2 
mo +m\ +ni2 

L 



\J k 2 (mo + m\)ai 



k 2 (mo -\-m\ +m2)fl2) 



(29) 
(30) 
(31) 



where m\ i are planetary masses, mo is the mass of the star, j = 
1,2 (we set the gravitational constant to k 2 ). To derive the octupole 
theory in terms of Jacobi reference frame, we start from writing 
down Eq. [2 with respect to Jacobi coordinates n 2 of two point 
masses, as a sum of two Keplerian terms and the perturbation: 



1 

2« 



Pi 



1 

2i"2 



m {) m 2 

INI 



-Hi 



where 



nf-l — 

P ert 



k 2 



m 1 mi 



r 2 -(i-Ki)n| 



-k m§m2 



1 

IN 



-Km 



Here, pi.2 are the conjugate Jacobi momenta, Ki = mi /(mo +m\ ), 
and the reduced masses are 



Mi 



mi mo 

(mo + mi ) ' 



("2 : 



m l{ m +mi) 
(mo +mi +m2) 



For details, see, e.g., (Mal hotnJll993h . Now, after expanding 
with respect to small Ki and retaining first order terms, we can show 
that the perturbation has the same form as Eq.[3] with the a ccuracy 
to the second order in the masses mi j2 /m dMalhotrdll993l) : 



-^in 



— k m\mi 



l|J- 2 — T! || ||r 2 || 

The s ame truncated Hamiltonian is analyzed in dLibert & Henrardl 
120051) who derived the secular Hamiltonian of the 12-th order in 
eccentricities by the "averaging with scissors" (i.e., by eliminating 
from the Fourier expansion of fC m all fast periodic terms depen- 
dent on /,■). These auth ors report that their secular t heory reproduces 
qualitatively results o f iMichtchenko & Malhotri] (2004) on the an- 
alytical way. 

The indirect part of the truncated Hamiltonian ?4« also aver- 
ages out to a constant, hence the rest of the averaging process is the 
same as in the case of f{ ptn written with respect to the Poincare el- 
ements. Nevertheless, the averaged !rt ml is missing terms of orders 
higher than two in the planetary masses. Indeed, with our method 
we derived the secular octupole Ha miltonian which ha s the same 
functional form as formulae (17) in jLee & Pealell2003l) . However, 
there are some differences in coefficients C2 and C3 [see their equa- 
tions (18) and (19)]. These coefficients in our expansion are the 
following: 



C 2 



C 3 



1 



(mo + mj 



16 (mo +m\ 
15 G 2 (m + 



|3, 



+ m 2) \mom\ 
m\^m\(mo — m\ 



3 ,3 r 3 D2 ' 



64 (mo + m\ +m2) 4 (momi) 5 L\G 2 



D 3 , 



(32) 



(33) 



where we extracted out two factors leading to the difference be- 
tween the respective formulae: 



D 2 



D 3 



1 



mo + mi 
mo 

(mp + mi) 2 
mo (mrj — mi ) 



O ^ 

\mo 

-l + O 



(34) 
(35) 
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These factors can be thought as equal to 1 in fcee & Pealell2003h . 
However, the theories are consistent within the assumed accuracy 
of the expansion and the relative magnitude of terms skipped from 
#j4t [of the order of 0(mi/mo)]. 

Actually, our averaging algorithm can be applied also to the 
full perturbing Hamiltonian, thanks to straightforward generaliza- 
tion for terms like the following: 

1 

||8iri -&2X2W ' 

where 81 , 82 are som e constants. In that instance, we obtain exactly 
the same C% 3 as in jLee & Pealell2003l) . Hence, as one would ex- 
pect, both approaches lead to fully equivalent results. 

This comparison also reveals that the averaging of the trun- 
cated Hamiltonian is in fact quite problematic because the accuracy 
of the secular expansion has nothing to do with the magnitude of 
the rejected terms. Already for Jupiter-mass planets, a contribution 
of these terms may be significant (see Sect. 2.4 for details). 

Moreover, a direct comparison of the theories would be more 
subtle. Although the functional forms of the perturbing Hamiltoni- 
ans are the same, they are expressed in terms of two different sets of 
canonical variables. Hence, the mean elements a,e,G5 have differ- 
ent meaning in these theories, i.e., the same physical configuration 
of the planets will be parameterized with quantitatively different 
values of the mean elements. 

2.4 Tests of the analytic secular theory 

To test the accuracy and relevance of the high order expansion of 
?4ec, we calculated the magnitude of subsequent terms relative to 
the free term. This expansion is computed up to the order of 24 for 
parameters (fij j,a,- ;) of extrasolar planetary systems taken from 
the Jean Schneider Encyclopedia of Extrasolar Planetfl The re- 
sults of this experiment are presented in Fig.Q] Each panel in this 
figure is labeled with the name of a relevant star and a number of 
putative planets it hosts (written in brackets). As we can see, for 
all examined systems, the sum of terms in J{ sec of the same order 
(note that we can have two and more planets in the system) decrease 
rapidly with the order of the expansion. For a few most separated 
systems with two planets (e.g., HD 217107, HD 190360), the high- 
est order terms are as low as ~ 10~ 37 in the relative magnitude. In 
the tested sample, the largest 24th-order terms are ~ 10~ 7 . Hence, 
the secular energy can be calculated with excellent accuracy. We 
note that sim ilar tests of the precision of the se cular theory were 
prese nted in dRodnguez & G allardo 2005) and I Liber t~& Henrardl 
l2005h . 

In the second test, we compare the outcome of our 
algorithm with the results o f semi-analytical a p proac h by 
Micht chenko & Malhotr~J d2004h; IMichtchenko et al.l ( 120061) [see 
also iMigaszewski & Gozdziewskil d2008al) for some technical as- 
pects] who applied the method to stud y the secular dynamics of 
two-pl anet system. In the algorithm of Michtchenko & M alhotral 
d2004l) the perturbing Hamiltonian is averaged out by means of 
the numerical integration. Hence, one avoids any expansion of the 
Hamiltonian and the results are formally exact (or very accurate, 
providing that precise enough quadratures are applied). 

To examine the accuracy of the analytic secular theory, we 
calculated the levels of the secular Hamiltonians for a number of 
extrasolar planetary systems which can be potentially regarded as 

1 http://exoplanets.eu 



non-resonant or well fitting assumptions of the secular theory. Be- 
cause the phase space of the secular problem is three-dimensional 
[i.e., (G\ , G2, AGS)], we choose the so called representative plane of 
the initial conditions to plot the secular energy levels. The defini- 
tion o f the representative plane follows IMichtchenko & Malhotral 
d2004l) . The secular Hamiltonian of a coplanar two-planet system 
depends only on ACS = G32 — tOj and eccentricities {e\,e%) cou- 
pled through the integral of the total angular momentum (or the so 
called Angular Momentum Deficit, AMD = G\ + G2). Effectively, 
the secular problem has one deg ree of freedom and is integrable. 
IMichtchenko & Malho tra (2004) have shown that all phase trajec- 
tories of the non-resonant s ystem pass throu gh a plane defined with 
AGJ = or AGS = 7t [see also dPauwelsfl9 83) for the qualitative anal- 
ysis of the secular two-planet problem]. The representative plane 
may be defined for fixed (X1.2 = a = a\ jai and fx\ 2 = ju = m\/ni2 
as follows: 

S = {e\ cosAG5 x e 2 ; ej e [0, 1), e2 £ [0, 1), AOS = U AGS = Jt}. 

This plane comprises of two (x = eicosAG5,y = £2)-half-planes 
with x ^ for AGS = 71 and with x > for AGS = 0. Simultane- 
ously, the derivatives of the secular Hamiltonian with respect to 
AGS are equal to zero for AGS = 0, 7t. It follows from the symme- 
try of interacting orbits with respect to both apsidal lines. Having 
the secular Hamiltonian in explicit analytic form, we can verify 
this property directly. Indeed, each term in the secular Hamiltonian 
depends on AGS only through cos(/AG5) with I e N, I > and it 
means that ^4 CC is even function of AGS. Hence, 3^/ scc /3ACi5 implies 
factors involving sin(ZAGS) and these terms vanish identically for 
AC! = 0,71. Formally, it is possible that 3.?4, C /3AG5 = also for 
AGS 7^ 0, 71, however, to find such solutions we should solve highly 
nonlinear equation involving (e\,e2) and trigonometric functions 
of AGS. 

Each pair of {e\,eq) for which dj{ xc /dG\ = corresponds to 
an equilibrium in the secular problem (simultaneously, they are the 
extrema of the secular Hamiltonian). These equilibria appear both 
in the negative half-plane of S , as mode II solutions (this mode is 
Lyapunov stable, and may be characterized with librations of an- 
gle AGS around 71 in the evolution of neighboring orbits), and in 
the positive half-plane of S as mode I solutions (Lyapunov sta- 
ble, with librations of angle AGS around of the nearby orbits). 
Further, in the regime of large eccentricities, in the positive half- 
plane a new, non-classic mode of motion may appear [it is the 
so called non-linear secular res onance, NSR from hereafter, see 
(Michtchenko & Malhotra 2004) for details]. These results are de- 
rived through the numerical (exact) approach, hence their reproduc- 
tion by the analytical theory provides an absolute test of its quality 
and accuracy. 

The results of the second experiment are presented in Figure|2] 
Each panel in this figure is for the representative plane of initial 
conditions computed and calculated for two-planet systems char- 
acterized with mass ratio fj and semi-major axes ratio a. These pa- 
rameters are written in the top-left corner at each relevant panel. To 
illustrate the significance of highest order terms in Eq. [22] we plot 
contour levels of (||#23|| + l-^^ID/H^bll, i- e -> the relative magni- 
tude of the sum of the last two terms with respect to the magnitude 
of the free term (in general, H n would stand for a sum of expan- 
sion terms over the number of planet pairs in the given multi-body 
configuration). Four particular contour levels of 1, 10 -8 , 10 -16 and 
10 -24 , respectively, are distinguished with thicker lines and labeled 
accordingly. 

In each panel, we also marked positions of stationary 
solutions (modes I, II, and NSR). Solutions represented by 
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Figure 1. Convergence of the expansion of the W-planet secular Hamiltonian derived in this paper. Each panel is for the magnitude of expansion terms 
in Eq. [22] divided by free term. The .v-axis is for the order of the expansion, the v-axis is for log| || W„||, where !H n denotes the magnitude of the sum of n-th 
order terms divided by the sum of the free terms. Parameters (/i; ,;,(X; ,) of ^4 ec are chosen for the known multi-planet extrasolar planetary systems. Their 
parent stars are marked in each panel together with number of planets written in brackets. 



thick curves are derived with the help of numerical averaging 
dMichtchenko & Malhotral [2004 ) : blue curves are for stable equi- 
libria, and thick green curves are for unstable solutions (the non- 
linear secular resonance). These curves represent practically exact 
(or very precise) solution to the problem. The thin, red lines mark 
the positions of equilibria calculated analytically, with the help of 
the 24th-order expansion of ^4 ec . 

The results can be summarized with a few interesting conclu- 
sions. Clearly, in the regions of the representative (e\ cosAtn,^)- 
plane at which (||#23|| + ||#24||)/||#b|| < 10~ 3 ,the precision of the 
analytical method is excellent. The secular theory predicts exactly 
positions of the equilibria in the negative half-plane of S (when 
AtO = 7r) no matter how large is ex. On contrary, the exact deriva- 
tion of the shape of the non-linear secular reson ance is a challeng- 
ing p roblem for the analytic approach (see also iHenrard & Libert! 
2005). The nonlinear resonance can be reproduced well by the 
analytic theory providing that (||#23|| + 1 1 ^24 1 1)/ 1 1^0 1 1 < 10~ 3 . 
This is also an empirical border of the convergence of the secu- 
lar expansion. We also found another empirical convergence con- 
dition that follows from the notion of the geometric series, i.e, 
11^24/^2211 < 1- This inequality is illustrated with triangular, or- 



ange colored regions labeled with "> 1". In these regions, the se- 
ries are divergent, hence the secular theory cannot reproduce the 
real dynamics. This may be interpreted as a clear limitation of the 
analytic theory. In fact, as we show below (Sect. 2.5), this problem 
appears rather due to imperfect algorithm of the expansion, which 
can be still improved. On the other hand, the results of this test 
provide an example that illustra tes excellent properties of the sem i- 
analytical approach invented by Mic htchenko & Malhotra (2004). 

2.5 An improved averaging algorithm 

A real source of the divergence of the secular series, Eq. [22] can 
be deduced after we draw in Fig.[2]the "anti-collision" line defined 
through 02(1 — e l) = a \ (1 + e l) ( see the right-half of the represen- 
tative plane). Clearly, the series diverge above this line and the po- 
sitions of the equilibria are strongly distorted. In this area, for some 
points or parts of the orbits, r,- [£,■(?)] > rj[fj(t)], while in the ex- 
pansion Eq.|22] r; < rj must be satisfied in the whole ranges of the 
anomalies. That may happen when the pericenter of the outer orbit 
is closer to the star than the apocenter of the inner orbit, no matter 
what is the relative orientation of their apsidal lines. Then the con- 
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Figure 2. A test of the secular theory of a coplanar two-planet system, derived in this paper. Each panel is for the representative energy plane, (e\ cos AH, £2)- 
In the right half-plane, AtD = 0, at the left half-plane, A03 = 7t. Gray region corresponds to crossing orbits and its boundary is defined through ai ( 1 — «2 ) = 
a\ (1 — e\ ), where (1) is for the inner orbit and (2) is for the outer orbit. Black thick line marks the "anti-collision" line defined with 02(1 — e-i) = ai(l + et) 
(see the text for more details). Black, thin lines are for contour levels of the last two terms of expansion Eq. 1221 i.e., the sum of absolute values of the 23-th and 
24-th order terms divided by the free term. A few contour levels (1, 10 -8 , 10 -16 and 10~ 24 , respectively) are marked with thicker curves. The area colored 
in orange determines the region where the expansion of 54 ec is divergent. Thick blue lines mark the positions of stable equilibria: mode II (with apsides 
anti-aligned, the left half-plane of the representative plan e), and mode I (with apsides alig ned, the right-half plane). These solutions are obtained numerically 
with the help of semi-analytical averaging algorithm in iMichtchenko & Malhotra 2004). Green curves are for the nonlinear- secular resonance (NSR). Red 
curves are for the corresponding libration modes calculated analytically with the help of the 24-order expansion of y/" S ec ■ Labels in the top-left corner at each 
panel are for the parameters of the expansion, a, = a\jai and p — mi/mi. 



dition of r, < r,; which required to write down Eq.Q~3] is violated. 
Obviously, it may be expressed by the equation of the anti-collision 
line, and in other words, through the requirement that the inner or- 
bit lies inside a circle of a radius equal to the pericenter distance of 
the outer orbit. Now it is also clear why for the apsides anti-aligned, 
we have always very good convergence of the secular series while 
in the right half-plane, the convergence region is generally strongly 
limited. Hence, the convergence limit of the expansion described 
in Sect. 2 may be simply interpreted through the conditions for 
the collision lines. We may also note that the problem persists in 
any secular theory that relies directly on Eq. [T3] One should be 



also aware that the conditions of crossing orbits do not depend on 
masses. If the masses are large, the dynamical collision curve ap- 
pears for much smaller eccentricity than the geometrical co llision 
line. (e.g jGozdziewski et alj20od : lMichtchenko et al]|2008h . 

A cure for the divergence problem may be a modified expan- 
sion of the term 1 /A/j, helping us to construct a secular theory that 
has no limits of the type + e;) < — ej). To derive the sec- 
ular expansion in Sect. 2, at first we factor rj in Eqs. ( lilt . (13). To 
have the series convergent for all positions of planets on their orbits, 
we propose to factor from the square root i I 1 [ 1 a term consisting of 
«2 multiplied by some scale factor, T| ^ 1 (instead of rj). Then we 



© 2008 RAS, MNRAS OOQ.mfKl 



Secular planetary problem 9 



can express the distance between planets ( and j as follows: 
A iJ = naj^l + ^,j(ri,rj), 



(36) 



where 



(37) 

The requirement of the convergent expansion of AfJ with re- 
spect to t^ij implies < 1. If = 1 the distance be- 
tween planets is equal to \/2r[aj. The convergence condition is 
fulfilled when maxA/j < \/2T]<Z/ for the given orbits. The maxi- 
mal distance between planets in coplanar orbits may be bounded 
by a, ( 1 + e; ) + a /(l + e /). Then the factor r| has the form of: 



n 



a;(l +ej) +aj(l +ej) 
\/2aj 



= -j=[aij(l + ei) + {l+ej)]. (38) 



For hierarchical systems with small eccentricity of the outer planet, 
T) ~ 1. For more compact systems with large eccentricities, T) ~ 
y2. For non-coplanar systems, T| may be even larger. In practice, 
this parameter makes it possible to control the convergence rate 
of the secular expansion. The convergence rate will be faster for 
large distances A but slower for smaller distances, moreover the 
condition of — 1 < j < 1 should be always fulfilled. 

The term 1 /Ay may be expanded with respect to f 



1 



1 

T\aj 



- (-i)'( 2 /-l)!! . 

2- 9/1 1 "»u 

/=! 



(39) 



To average out the above formulae, we express positions of both 
planets in a given pair with respect to the eccentric anomaly. Next, 
we change the integration variables similarly as in Sect. 2, i.e., 
dMk = I^iE^e^jdE^. We also should express r,- ■ tj through AB3 
(again, fixing the reference frame with the apsidal line of the inner 
orbit). After expressing terms r,-, ri , r, ■ ti through eccentric anoma- 
lies, the function £ = £12 has the following explicit form (to shorten 
the notation, let us fix i = 1, j = 2 for a given pair of planets): : 

C= -S T |e + eicos£i+e2Cos£2 + 63sin£'i+ (40) 
r\ z L 

+64 sini?2 + 85 cos£i 2 + 9g cos£2 2 + 67 cos£i cos £2 + 
+6g s'mE\ sin£2 + 69 cos£"i sin£2 + 6io cos ^2 sin£i , 



where E\ ,£2 are eccentric anomalies of the inner and outer planets 
respectively, and coefficients 0/,/ 0, read as follows: 

,2 



6o 
6| 

62 

63 



: a~ — 2ae 1 e2COsA05+ 1 — T| 
2a<?2 cosAG5-2a 2 ei, 
: 2ae! cosACJ — 2e2, 

■■ -2ae 



e\ sin ACS, 



= 2aei\l\ 



e 4 

e 5 
e 6 

e 7 = — 2acosAro, 



e\ sinACB, 



2 2 
■el 



(41) 



-2a\/i — e\J\ — e\ cos ACS, 



e 9 

810 = 



= — 2fXi / 1 — ei sin ACS, 



:2a\/l 



e\ sinAtO. 



Here, a = ttj^, AtD = AG5i2 and ei,«2 are the eccentricities of the 
inner and outer planet, respectively. 

After the double averaging of X, over the mean anomalies we 
obtain: 



J-j [rj + r) 2r, • r,- - n^] - ^2 [ A tj ■ < ? >= 



1 

2rj2 



3e 2 + 2j a 2 -9aeie2CosAC3 + 3e2 -2r) 2 +2 



The averaging of the square of £ over the mean anomalies brings 
the following formulae: 



1 



a 4 (8 + 40e 2 + 15et 



8r| 4 

a 3 ^-30eie 2 (3e 2 +4)cosA05) + (42) 
a 2 [(2-ri 2 )(16 + 24e 2 )+e 2 e|(72+100cos2A0i5)+48e^ + 



-a( -6eie 2 (20- 12r| 2 + 15e?;)cosA0J) + 



-8(ri 2 - l) 2 +40e2-24e^r| 2 + 154 



These preliminary calculations show that the new algorithm leads 
to more complex expansion of the secular Hamiltonian than the 
simple approach in Sect. 2.2 which, as we have demonstrated, is 
limited in some cases. Moreover, we found this improvement af- 
ter submitting the manuscript, hence the new expansion and a de- 
tailed study of its properties would make the paper very lengthy. We 
are going to present the improved algorithm and the results of its 
tests in a new work devoted to the analytic theory of non-coplanar 
model of W-planets. A generalization of Eq.[38]for that case seem 
straightforward, because we should only calculate r, • rj with the 
help of appropriate rotation matrix parameterized through Euler an- 
gles, i.e., the Keplerian elements (ii t ij,(Oi, G>j,ili,£lj). Hence, only 
the coefficients 0/ will be modified. 



3 SECULAR DYNAMICS OF THREE-PLANET SYSTEM 

The two-planet secular Hamiltonian in Sect. 2 can be easily adapted 
to construct the secular theory for A'-planet system. At present, 
a few candidates of such configurations are a lready discovered , 
including four planet systems, e.g., y Arae dJones et al ] 120021: 



iButler et al.l2006al : lGozdziewski et al.l2007l : |Pepe et alj2007l). five 
or ev en six planet configuration around 55 Cnc dFischer et al.l 
2008). Here, as the simplest and most natural generalization of the 
two-planet model, we consider the secular theory of three-planet 
configuration which is far from MMRs and collision zones. 

The three-planet model is described by Hamiltonian in 
Eqs. |Q3, $2% and <[3j, respectively, where N = 3. Because the nodal 
longitudes are undefined in the coplanar system, and the secu- 
lar dynamics depends on the relative positions of the mean or- 
bits, we can eliminate the nodal longitudes from the problem. Let 
indices i = 1,2,3 enumerate the planets. Their semi-major axes 
are a\ < a-i < 03, respectively. After averaging 91 over the mean 
anomalies, the secular Hamiltonian #" sec does not depend on /, any- 
more. Therefore, conjugate momenta L, (hence, the semi-major 
axes) are constants of motion. Because the secular system does not 
depend on particular longitudes of nodes, the respective degree of 
freedom is also irrelevant for the secular dynamics. 

Hence, the secular system can be described with the following 
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set of canonical elements: 



gl = 


-BJl, 


Gi, 




82 = 


-052, 


G 2 , 


(43) 


83 = 


-ra 3 , 


G 3 . 





We can eliminate one more degree of freedom with the help of 
the angular momentum integral, which can be also expressed with 
AMD = Gi + Gi + G3. For that purpose, we perform the following 
canonical transformation: 



CT l =81 -83 = ®3 —GIi = Ado, Gi, 
^2 = g2-g3 = ^3 -ra 2 = Am 2 . 3 , G 2 , 
0-3 =g 3 = -C33, AM£» = Gi+G 2 + G3, 



(44) 



introducing new canonical angles 01,02,03. These angles can be 
interpreted as two-planet ACS defined for each pair of planets in 
the three-planet system. The secular Hamiltonian can be expressed 
through Oi and o 2 explicitly, hence 03 is cyclic and AMD is con- 
stant of motion. Actually, we reduced the secular system of three 
planets to two degrees of freedom, with the secular energy and 
AMD as free parameters. 

3.1 Representative planes of initial conditions 

Now, we have a similar problem as in the case of two-planet con- 
figuration. We want to illustrate the dynamical properties of the 
system in possibly global manner. To characterize its dynamical 
states, we follow the general idea of the representative plane of 
initial conditions. Here, we focus on the equilibria in the secular 
problem. Fixing the integral of AMD as a parameter of the system, 
the dynamics may be represented in four-dimensional phase space 
of (Gi , G2,Oi ,02), or (ei ,e2,^i ,CT 2 ). The representative plane will 
be chosen according with: 



= 0, 



0. 



(45) 



3cji 3c*2 

Then any pair of points (e^e®) that belongs to the representative 
plane, and the following equations are satisfied: 



= 0. 



:0, 



(46) 



3Gi 9G2 
defines an equilibrium of the secular problem. Note that the last 
conditions implies also 3^/ scc /3G3 = because G3 = G3(Gi,G2) 
is a function parameterized by the total angular momentum (or 
AMD). The explicit and simple transformation between the eccen- 
tricity and the element G makes it possible to solve the above con- 
ditions in terms of (e\ , e 2 ). 

Actually, the most obvious definition of the representative 
plane is a generalization of that plane constructed in the two-planet 
problem. For fixed a 1 , a 2 , 03 , in 1 , m 2 , 1113 and AMD, as a parameter, 
the symmetric representative plane is the set of points such that 

5 = {e\ cos 0i x e2coso 2 ; e \ £ [0> 1), e 2 £ [0, 1), CJ12 = OU71}. 

This plane comprises of four quaters. The signs of ej 2 of the coor- 
dinated axes tell us on the respective values of the secular angles. 
In that case, the condition in Eq. [45] is fulfilled thanks to apsidal 
symmetries of the problem. It is also obvious by recalling that Jf sec 
is even function of o"i. 2 . The derivatives of J4 CC over o"i. 2 depend on 
factors involving sin(/ai.2), / £ N, / > that vanish for Oi 2 = 0,71. 
Alternatively, the representative plane may be also defined through 
sine 12 = 0. 

Moreover the condition for the representative plane, defined 



through the vanishing derivatives over secular angles, may be also 
fulfilled for sino"i. 2 7^ 0. Let us start with the octupole (third-order) 
approximation of the secular Hamiltonian. Using Eqs. {22}, ( 124b 
and d25K we can write the secular Hamiltonian in the following 
short form of: 

#see =Yl,2COsAn3i,2 + Yi,3COsAG5i.3+Y2.3COsAG3 2 ,3+Y4- (47) 
Using the canonical angles defined with Eq.1441 

#sec = Yi,2Cos(oi -<t 2 )+Yi,3Cos<Ti +Y2.3COSO2+Y4, (48) 

where Yl.2, Yl.3> Y2.3, Y4 Me functions of e\,e 2 ,ej. Hence, the non- 
symmetric representative plane is defined through the following 
conditions, the same as Eq.[45] in the explicit form: 

-Yi,2sin(ai -02) -Yusinrji =0, 

Yi,2sin(fJi — 02) -Y2,3 sm °2 =0. (49) 

Obviously, conditions in Eq. [49] are satisfied not only when 
sin<Ti,2 = 0. We have four other solutions, satisfying Eq.[49] i.e.: 



o"i = iarccos 



Oi = ± arccos 



1 



, | 72,3 Y2,3 

2 1,2 I Y2,3 Y?2 Y13 



1 



1 

Yl,2 



Yl,2 71,2 

To Y23 



(50) 



These solutions describe the non-symmetric representative planes, 
with respect to the octupole theory. This definition is exact up to 
the third order in (X; /. Angles Oi , o"2 satisfying condition in Eq.1451 
i.e., solutions to Eq. [49] can be found consistent with higher or- 
der expansions. However, these equations are very complex and, in 
practice, we would have to solve them numerically (for instance, 
with the Newton-Raphson algorithm initiated with starting condi- 
tions derived from the octupole theory). Hence, in general, the rep- 
resentation of the energy levels with the help of the non-symmetric 
representative planes is much more difficult than in the symmet- 
ric case and is not unique (it has some analogy to the Poincare 
cross-section). For instance, we can define the representative plane 
for higher order expansions of ^ scc . In this paper, we focus on the 
symmetric representation only. 

3.2 Energy levels for three-planet secular model 

In this section, to show some applications of the secular theory, 
we investigate qualitative dynamics of a few three-planet extraso- 
lar systems. To characterize these systems, we calculated energy 
levels in the symmetric representative plane. We also try to find 
equilibria in the secular model of each examined system. The re- 
sults are illustrated in Figure [3] We selected four three-planet con- 
figurations. Their orbital elements are taken from the Extrasolar 
Planets Encyclopedia of Jean Schneider, following the most re- 
cent determinations of the orbital solutions. E ach column in Fig. [3 
is for one p articular system, i.e., for rj And (Butler e t alj|20061 



HD 74156 (Beanetal 



HD 69830 (Udr vetal 



2008), Gliese 581 dLovis et alj 2006) and 
2007), respectively. Different approxima- 
tions to the secular theory are illustrated in rows. The top row pan- 
els are for the energy levels calculated with the octupole theory, 
panels in the middle row are derived from the expansion of J{ xc 
of the 24-order, and panels in the bottom row illustrate the energy 
levels computed with the numerical algorithm. The later case may 
be regarded as the exact solution to the problem thanks to adaptive, 
high order Gauss-Legendre quadratures which we used to compute 
the double integral over ^{ pat . In each case, we fixed AMD consis- 
tent with the nominal parameters of the examined systems. 



© 2008 RAS, MNRAS 000,[TJ{T4] 



Secular planetary problem 1 1 




Gl5t 


1 

4. 


oncuyticol 


\ 




1 



-0.4 -0.2 0.0 0.2 0.4 





-1.0 -0.5 0.0 0.5 



1.0 -0.5 0.0 0.5 



-0.4 -0.2 0.0 0.2 0.4 



-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 



Figure 3. The secular energy levels on the symmetric representative plane for selected three-planet systems. Map coordinates are (x = e\ coso~i , y = £2 cos (5t), 
where the secular angles o~i,o~2 are or %. Black thin lines are for energy levels obtained with different methods: panels in the top row are for the octupole 
theory, panels in the middle row are for the 24th-order expansion derived in this paper, and panels in the bottom row are for the semi-analytical averag- 
ing. In gray areas, £3 < 0, hence the motions are not permitted The light-gray areas are for the regions of collisions between the inner and the middle 
planet. In orange regions, the secular expansion diverges. The thick, red straight lines mark the anti-collision lines between the planets and indicate the 
true border of the convergence of the secular expansion. Orbital parameters are taken from Jean Schneider Encyclopedia, and are given in terms of tuples 
p, = (mo[M ],mi[mj],m2[mj],m 3 [mj],fl 1 [AU],O2[AU],<33[AU],e 1 ,e2,e3) as follows: p upsA nd = (1.27,0.69,1.98,3.95,0.059,0.83,2.51,0.029,0.254,0.242), 
Phd 74156 = (1-24, 1.88,0.396,8.03,0.294, 1.01,3.85,0.64,0.25,0.43), p G i58i =(0.31,0.0492,0.0158,0.0243,0.041,0.073,0.25,0.02,0.16,0.2), Phd 69830 = 
(0.86,0.033,0.038,0.058,0.0785,0.186,0.63,0.1,0.13,0.07). Each panel is labeled with AMD (expressed in standard units) calculated for the nominal con- 
figuration. The secular energy levels for the nominal systems are marked with green curves. 



3.2.1 V Andromedae 



First, we are looking at the exact (numerical) phase plots. In the 
case of D And, we found two types of equilibria. The first one, 
marked with I, is the global minimum of the secular Hamiltonian. 
Hence, according to the Lyapunov theorem, this equilibrium is sta- 
ble (because the Hamiltonian can be regarded as the positive defi- 
nite Lyapunov function). The equilibrium marked with II is a saddle 
point, and is stable in the linear approximation. It can be verified 
by solving the eigenproblem of the linearized equations of motion 
in the neighborhood of the equilibrium. Now, we can compare the 
outcomes of the analytic theories with the results of exact, numeri- 
cal algorithm. Apparently, the high-order analytical theory is fully 
compatible with the numerical theory. The largest deviations be- 
tween the secular energies are of the order of 10~ 9 . This accuracy 
is preserved even for eccentricities e\ close to 1. On contrary, the 
octupole theory provides only a crude representation of the phase 
space. The phase plot constructed with the help of this theory also 
reveals an equilibrium of type I, however, at place of equilibrium II, 



qualitatively different energy levels appear (see the top-left panel in 
Fig.[3]with three "false" equilibria). 

To locate the "real" system in the energy plot, we mark the 
level of the secular energy computed for the nominal parameters of 
the system with green, thick curve. It provides only a crude imag- 
ination where the system is located; one should be aware that we 
are looking at the representative plane (hence C\ 2 are fixed at spe- 
cific values), and we do not take into account the parameter errors. 
Still, the plot tell us that while variability of is limited, e\ may 
be varied in all permitted range of eccentricity. 

3.2.2 HD 74156 

In the phase space of the HD 74156 system, we discover only one 
equilibrium (labeled with III) in the regime of small eccentricities. 
It is related to the global maximum of ^ scc , and it means that this 
solution is Lyapunov stable. The AMD of the nominal system per- 
mit eccentricities to reach large values, hence they enter the regions 
in which the secular Hamiltonian expansion diverges (see the ex- 
planation in Sect. 2.5). These regions are marked in orange color. 
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The region of permitted motions is also bordered by two collision 
lines of orbits, defined implicitly through az(l — ei) = ai(l ± e\) 
and they are marked with gray, thick lines. In this case the view 
of the phase plot varies with the order of expansion (or the applied 
algorithm). The high-order analytic theory reconstructs the phase 
plot for «2 ~ [0, 0.6) (white zone) in almost whole permitted range, 
nevertheless, in the region of divergent expansion the phase plot 
is wrong. In the case of octupole theory, we obtain only a crude 
approximation of the structure of the phase space, and again the 
theory introduces artifacts (two saddle points and an extremum). 

We also plot the energy levels of the nominal system (in the 
same manner as we did for l) And). Its parameters would evolve 
along this level relatively distant from the equilibrium close to the 
origin. 

3.2.3 Gliese 581 and HD 69830 

The phase space of systems Gliese 581 and HD 69830 are quite 
similar. The region of permitted motions is limited to relatively 
small eccentricities. In both cases, we have only one equilibrium 
close to the origin that is related to the global maximum of the sec- 
ular Hamiltonian, and therefore they are Lyapunov stable. For the 
Gliese 581 planetary system, the accuracy of the 24-order secular 
expansion is not very good; at the borders of permitted motion, 
this accuracy is at a level of 10~ 3 only. For the third-order the- 
ory, this accuracy is even worse, ~ 10~ 2 . In the case of HD 69830, 
the relative accuracy of the high-order expansion is not worse than 
2.5 x 10~ 9 . 

3.3 Secular dynamics of HD 37124 

As a particular system to study, we choose th e three-planet sys - 
tem of HD 37124. It has been discovered by IVogt et al.l J2005h . 
Remarkably, the most recent best fit solutions to the observations 
are consistent with configurations involving sub-Jupiter compan- 
ions in orbits with moderate eccentricities. The eccentricity of the 
outermost companion is not well constr ained, nevertheless exten- 
sive dynamical analysis of the RV data in ( Gozdziewski et al.l2008h 
make it possible to locate this planet in a region between 8:3 and 
11:4 MMRs with the middle companion. In that case, the orbital 
parameters can be regarded as well fitting the assumptions of the 
secular theory. 

Figure [4] is for the energy levels in the symmetric represen- 
tative plane computed and calculated for three slightly different 
orbital configurations related to possible orbital best-fits (panels 
in each column are for one orbital fit). Osculating elements of 
these configurations are quoted in the caption to this figure. The 
first, kinematic solution to the three- Keple rian model of th e RV, is 
taken from the original dis covery paper jVogt et alj [2005). Other 
two best-fits are from lGozdziewski et al.l J2008h . The energy levels 
are computed for fixed AMD calculated for each particular initial 
condition. 

In the top row of Fig. [4] we show energy levels calculated from 
the octupole theory, plots in middle row are derived from the ex- 
pansion of of the 24-order, and the bottom row illustrates the 
results derived from the numerical algorithm. In all cases, the ana- 
lytical high-order theory is in excellent agreement with the numeri- 
cal, exact theory. We checked that the magnitude of largest, relative 
deviations between the analytic and numerical results are of the or- 
der of 2.5 x 10~ 6 . On the other hand, the octupole theory gives 
relatively precise insight into the structure of the phase space. All 
qualitative features of the energy plane are reproduced quite well. 



The HD 37124 seems to be the most interesting example of the 
secular dynamics in the real system found in this paper. The energy 
planes reveal unusual dynamical structures related to the equilibria 
in the secular system. We know already that they can appear as ex- 
trema as well as saddle points in the representative plane. For the 
same AMD, we can have three types of stationary solutions. Two of 
them are characterized by extrema of J{ sec in the four-dimensional 
phase space: the first one is the maximum which appears in the 
quarter with (J1.2 = 0, and there is the global minimum of J4ec in 
the quarter with G\ = and 02 = 31. We also found saddle points of 

in the quarter with G\ = 7t, 02 = 0, Jt. Stationary points marked 
with I and II in Fig.|4]are related to extrema of hence they are 
stable. By examining the eigenvalues of the linearized equations of 
motions in the neighborhood of the saddle point (equilibrium III), 
we checked out that it is linearly stable. It can be localized in the 
half-plane of o"2 = or C2 = 1, depending on selected orbital pa- 
rameters. 

All these solutions appear in the range of moderate eccentric- 
ities, and in fact can be located close to the actual positions of the 
best fit solutions. The structure of the energy plane is also robust 
with respect to small changes of the orbital parameters. In each 
panel, similarly to the previous systems, we mark the energy level 
of the respective nominal configuration with the green thick curve. 
Curiously, depending on the chosen fit, the nominal system can 
evolve in the quarter of the representative plane characterized by 
librations of <5\ 2 around (the top-right quarter), or librations of 
Oi around and 02 around Jt (the bottom-right quarter), as well as 
C\ around 71 while 02 can be librating around or n (the top-left or 
the bottom-left quarter). It means, that the apses of two innermost 
companions can be all aligned with the apsidal line of the outer- 
most planet, they can be also anti-aligned or the apsidal directions 
can be mixed. 

The presence of these stationary solutions can be interpreted 
in terms of the two-planet theory. We recall that for the case of 
two planets, mode I (with apsides aligned) corresponds to the max- 
imum of the secular energy, while mode II (apsides anti-aligned) 
corresponds to the minimum of ?{ scc . In the case of three planets, 
we add the secular energies of three pairs of interacting planets. 
Hence, in a region of the representative plane where the maxima 
of J/sec of these three planetary pairs can roughly coincide, we can 
obtain the maximum of the total energy; by adding ^4 CC in the re- 
gion where the particular minima are close enough in the parameter 
space, we can obtain the global minimum of the energy, and in the 
case of superimposed minimum and other maximum we can obtain 
the saddle of the total energy. Geometrically, the equilibria can be 
interpreted as combinations of the secular modes known from the 
theory of nonresonant two-planet system. For instance, the maxi- 
mum of 54ec can be related to triple mode I (i.e., the neighboring 
solutions are characterized with librations AtD around for all pairs 
of planets), and the saddle point of J{ XQ is obtained for a super- 
position of mode I for some pair(s) and of mode II for the other 
pair(s). It is not clear for us yet, what would mean a combination 
with the NSR mode, in the regime of large eccentricities (i.e., in the 
region of nonlinear-secular resonance). Likely, it could be related 
to sophisticated secular dynamics. 



4 CONCLUSIONS 

The number of multi-planet extrasolar systems constantly grows. 
The Doppler spectroscopy remains the most effective detection 
technique. Unfortunately, the measurements of RV are in some 
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Figure 4. The secular energy levels on the symmetric representative plane for three-planet system HD37124 (see the text for more details). The map coordinates 
are (x=e\ cosCi ) and (y = ei coso?), where the secular angles are fixed at (the top half-plane/the right half-plane) or % (the bottom half-plane/the left half- 
plane). The boundary between white and grey regions is for ej = (i.e., the motion is permitted only in the white areas). Black solid lines are for the secular en- 
ergy levels obtained with the help of the octupole theory (panels in the top row), by the expansion of the 24-th order (panels in the middle row) and with the semi- 
numerical averaging (the bottom row of panels). Orbital parameters are taken from the discovery paper (Vogt et al., 2005) (panels in the left column marked 
with 1), and from Gozdziewski et al. (2007) (panels in the middle and in the right column, marked with 2 and 3, respectively). These parameters, in terms of tu- 
ples pj = (wo[M0],wi [mj],m2[mj],m 3 [mj],ai[AU],a2[AU],fl3[AU],ei,e2,e3) are the following: p] = (0.91,0.61,0.6,0.683,0.53,1.64,3.19,0.055,0.14,0.2), 
p 2 = (0.78,0.624,0.606,0.581,0.519,1.632,3.212,0.037,0.003,0.048), p 3 = (0.78,0.650,0.584,0.567,0.519,1.668,2.740,0.091,0.040,0.132). Stationary 
solutions are labeled with I, II, and III. Green curves are for the secular energy level of the respective nominal initial condition. 



sense degenerate because due to symmetry of the Doppler signal, 
we usually cannot determine the true inclination of planetary orbits. 
Other parameters are usually determined with large uncertainties. 
Hence, to characterize the dynamics of such systems we cannot rely 
only on single initial conditions and effective, qualitative methods 
of dynamical analysis are very desirable. 

In this paper we consider the secular theory of a copla- 
nar iV-planet system which is far from MMRs and orbital col- 
lision zones. In this case, the high-frequency interactions can 
be averaged out and we obtain greatly simplified picture of the 
long-term behavior of the system. This idea leads to the classic 
Laplace-Lagrang e theory and its modern generalizat ions like the 
octupole theory dFord et aljbood : iLee & Pealell2003l). high-order 
expansion of the secular perturbation drlenrard & Libertl 1 20051 : 



iRodrfguez & Gallarddl2005h or the semi-n umerical averaging in- 
vented bv lMichtchenko & Malhotr3 d2004h . 

Our work can be considered as a generalization of the octupole 
theory for hierarchical triple systems characterized with large ratio 
of semi-major axes. We have shown that in this case the pertur- 
bation can be averaged out over mean longitudes with very basic 
change of integration variables that makes it possible to express 
the integrand function as a polynomial of trigonometric functions, 
without any need of relatively complex Fourier expansions. To the 
best of our knowledge, such method has been not applied in the lit- 
erature. However , during the final preparation of the manuscript we 
found a book of [Valtonen & Karttunen] d2006j) . who use a similar 
idea to construct the octupole theory of hierarchical triple stellar 
systems. 

Basically, the secular Hamiltonian is expressed through the 
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power series with respect to the ratios of semi-major axes, without 
an explicit limitation on the eccentricities. These series can be con- 
tinued to practically any order. However, if we apply the averaging 
algorithm presented in this work, the convergence region of these 
series is usually limited. To avoid this problem, in this paper we also 
propose a further improvement of this method and to generalize the 
expansion to the spatial problem. Our theory significantly improves 
the octupole theory of two-planet systems. We have shown that it 
can be generalized to any A'-planet system fulfilling the assumption 
of the averaging theorem. The simple "trick" of choosing the inte- 
gration variables can be applied not only for purely gravitational 
point-to-point interactions but also in other models in which the 
mutual interactions can be expressed in powers of mutual distance 
between objects in the system. For instance, now we work on apply- 
ing the averag ing algorithm to relativistic and quad rupole moment 
perturbations dMigaszewski & Gozdzie wski 2008b). Its generaliza- 
tion to the 3D problem (in particular, for two-planet system) is also 
straightforward. Then it can be applied to the study of secular dy- 
namics in hierarchical triple-star systems or star-planet configura- 
tions fulfilling assumptions of the secular theory. 

In this work, the secular theory is used to investigate station- 
ary solutions in the three-planet systems that are relatively frequent 
in the known sample of extrasolar planets. We found that the libra- 
tion modes known in two-planet configurations can be generalized 
for the multi-planet model. Still, our study of particular systems 
is quite preliminary and new, yet unknown stationary solutions are 
expected to exist in this problem. 
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